home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / bessel.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  35.9 KB  |  1,179 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8.  
  9. (in-package "MAXIMA")
  10.  
  11. ;; When non-NIL, the Bessel functions of half-integral order are
  12. ;; expanded in terms of elementary functions.
  13. (defmvar $besselexpand nil)
  14.  
  15. ;; Temporarily we establish an array convention for conversion
  16. ;; of this file to new type arrays.
  17.  
  18. (eval-when (compile eval)
  19.  
  20. ;; It is more efficient to use the value cell, and we can probably
  21. ;; do this everywhere, but for now just use it in this file.
  22.        
  23. (defmacro nsymbol-array (x) `(symbol-value ,x))
  24. ;(defmacro nsymbol-array (x) `(get ,x 'array))
  25.  
  26. (defmacro narray (x typ &rest dims) typ
  27.   `(setf (nsymbol-array ',x)
  28.      (make-array
  29.        ,(if (cdr dims) `(mapcar '1+ (list ,@ dims))
  30.           `(1+ ,(car dims))))))
  31. )
  32.  
  33. (declare-top(flonum (j[0]-bessel flonum) (j[1]-bessel flonum)
  34.          (j[n]-bessel flonum fixnum) (i[0]-bessel flonum)
  35.          (i[1]-bessel flonum) (i[n]-bessel flonum fixnum)
  36.          (g[0]-bessel flonum) (g[1]-bessel flonum)
  37.          (g[n]-bessel flonum fixnum))
  38.      (flonum x z y xa sx0 sq co si q p)
  39.      (special $jarray $iarray $garray)
  40.      (array* (flonum j-bessel-array 1. i-bessel-array 1.
  41.              g-bessel-array 1.))
  42.      (array* (flonum $jarray 1. $iarray 1. $garray 1.))
  43.      (*fexpr $array)) 
  44.  
  45. #-(or cl NIL)
  46. (and (not (get '*f 'subr)) 
  47.      (mapc #'(lambda (x) (putprop x '(arith fasl dsk liblsp) 'autoload))
  48.        '(*f //f _f +f -f)))
  49.  
  50. #-NIL
  51. (declare-top(flonum (*f flonum flonum) (//f flonum flonum) 
  52.          (_f flonum fixnum) (+f flonum flonum) (-f flonum flonum))
  53.      (*expr *f //f _f +f -f))
  54.  
  55. #+(or cl NIL)
  56. (eval-when (eval compile)
  57.   (defmacro *f (a b) `(*$ ,a ,b))
  58.   (defmacro //f (a b) `(//$ ,a ,b))
  59.   (defmacro +f (a b) `(+$ ,a ,b))
  60.   (defmacro -f (a b) `(-$ ,a ,b))
  61.   ;_f isn't used here.  That would be scale-float, no open-code version.
  62.   )
  63.  
  64.  
  65. ;;
  66. ;; Bessel function of the first kind of order 0.
  67. ;;
  68. ;; One definition is
  69. ;;
  70. ;;         INF
  71. ;;         ====       k  2 k
  72. ;;         \     (- 1)  z
  73. ;;          >    -----------
  74. ;;         /       2 k   2
  75. ;;         ====   2    k!
  76. ;;         k = 0
  77. ;;
  78. ;; We only support computing this for real z.
  79. ;;
  80. (defun j[0]-bessel (x) 
  81.    (slatec:dbesj0 (float x 1d0)))
  82.  
  83. (defun $j0 ($x)
  84.   (cond ((numberp $x)
  85.      (j[0]-bessel (float $x)))
  86.     (t (list '($j0 simp) $x))))
  87.  
  88.  
  89. ;; Bessel function of the first kind of order 1.
  90. ;;
  91. ;; One definition is
  92. ;;
  93. ;;      INF
  94. ;;      ====   - 2 k - 1      k  2 k + 1
  95. ;;      \     2          (- 1)  z
  96. ;;       >    --------------------------
  97. ;;      /            k! (k + 1)!
  98. ;;      ====
  99. ;;      k = 0
  100.  
  101. (defun j[1]-bessel (x) 
  102.    (slatec:dbesj1 (float x 1d0)))
  103.  
  104. (defun $j1 ($x)
  105.   (cond ((numberp $x)
  106.      (j[1]-bessel (float $x)))
  107.     (t (list '($j1 simp) $x))))
  108.  
  109. ;; Bessel function of the first kind of order n
  110. ;;
  111. ;; The order n must be a non-negative real.
  112. (defun $jn ($x $n)
  113.   (cond ((and (numberp $x) (numberp $n) (>= $n 0))
  114.      (multiple-value-bind (n alpha)
  115.          (floor (float $n))
  116.        (let ((jvals (make-array (1+ n) :element-type 'double-float)))
  117.          (slatec:dbesj (float $x) alpha (1+ n) jvals 0)
  118.          (narray $jarray $float n)
  119.          (fillarray (nsymbol-array '$jarray) jvals)
  120.          (aref jvals n))))
  121.     (t (list '($jn simp) $x $n))))
  122.  
  123.  
  124. ;; Modified Bessel function of the first kind of order 0.  This is
  125. ;; related to J[0] via
  126. ;;
  127. ;; I[0](z) = J[0](z*exp(%pi*%i/2))
  128. ;;
  129. ;; and
  130. ;;
  131. ;;        INF
  132. ;;        ====         2 k
  133. ;;        \           z
  134. ;;         >    ----------------
  135. ;;        /         2 k   2 
  136. ;;        ====     2    k!
  137. ;;        k = 0
  138.  
  139. (defun i[0]-bessel (x)
  140.    (slatec:dbesi0 (float x 1d0)))
  141.  
  142. (defun $i0 ($x)
  143.   (cond ((numberp $x)
  144.      (i[0]-bessel (float $x)))
  145.     (t (list '($i0 simp) $x))))
  146.  
  147. ;; Modified Bessel function of the first kind of order 1.  This is
  148. ;; related to J[1] via
  149. ;;
  150. ;; I[1](z) = exp(-%pi*%I/2)*J[0](z*exp(%pi*%i/2))
  151. ;;
  152. ;; and
  153. ;;
  154. ;;       INF
  155. ;;       ====         2 k
  156. ;;       \           z
  157. ;;        >    ----------------
  158. ;;       /      2 k
  159. ;;       ====  2    k! (k + 1)!
  160. ;;       k = 0
  161.  
  162. (defun i[1]-bessel (x)
  163.   (slatec:dbesi1 (float x 1d0)))
  164.  
  165. (defun $i1 ($x)
  166.   (cond ((numberp $x) (i[1]-bessel (float $x)))
  167.     (t (list '($i1 simp) $x))))
  168.  
  169. ;; Modified Bessel function of the first kind of order n, where n is a
  170. ;; non-negative real.
  171. (defun $in ($x $n)
  172.   (cond ((and (numberp $x) (numberp $n) (>= $n 0))
  173.      (multiple-value-bind (n alpha)
  174.          (floor (float $n))
  175.        (let ((jvals (make-array (1+ n) :element-type 'double-float)))
  176.          (slatec:dbesi (float $x) alpha 1 (1+ n) jvals 0)
  177.          (narray $iarray $float n)
  178.          (fillarray (nsymbol-array '$iarray) jvals)
  179.          (aref jvals n))))
  180.     (t (list '($in simp) $x $n))))
  181.  
  182. (defun $bessel_i (arg order)
  183.   (if (and (numberp order)
  184.        (numberp ($realpart arg))
  185.        (numberp ($imagpart arg)))
  186.       ($in (complex ($realpart arg) ($imagpart arg)) order)
  187.       (subfunmakes '$bessel_i (ncons order) (ncons arg))))
  188.  
  189. (defun bessel-i (arg order)
  190.   (cond ((zerop (imagpart arg))
  191.      ;; We have numeric args and the first arg is purely
  192.      ;; real. Call the real-valued Bessel function.  We call i0
  193.      ;; and i1 instead of jn, if possible.
  194.      (let ((arg (realpart arg)))
  195.        (cond ((zerop order)
  196.           (slatec:dbesi0 arg))
  197.          ((= order 1)
  198.           (slatec:dbesi1 arg))
  199.          (t
  200.           (multiple-value-bind (n alpha)
  201.               (floor (float order))
  202.             (let ((jvals (make-array (1+ n) :element-type 'double-float)))
  203.               (slatec:dbesi (float (realpart arg)) alpha 1 (1+ n) jvals 0)
  204.               (narray $besselarray $float n)
  205.               (fillarray (nsymbol-array '$besselarray) jvals)
  206.               (aref jvals n)))))))
  207.     (t
  208.      ;; The first arg is complex.  Use the complex-valued Bessel
  209.      ;; function
  210.      (multiple-value-bind (n alpha)
  211.          (floor (float order))
  212.        (let ((cyr (make-array (1+ n) :element-type 'double-float))
  213.          (cyi (make-array (1+ n) :element-type 'double-float)))
  214.          (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
  215.                     v-cyr v-cyi v-nz v-ierr)
  216.          (slatec::zbesi (float (realpart arg))
  217.                 (float (imagpart arg))
  218.                 alpha
  219.                 1
  220.                 (1+ n)
  221.                 cyr
  222.                 cyi
  223.                 0
  224.                 0)
  225.            (declare (ignore v-zr v-zi v-fnu v-kode v-n
  226.                 v-cyr v-cyi))
  227.  
  228.            ;; We should check for errors here based on the
  229.            ;; value of v-ierr.
  230.            (when (plusp v-ierr)
  231.          (format t "zbesi ierr = ~A~%" v-ierr))
  232.            (narray $besselarray $complete (1+ n))
  233.            (dotimes (k (1+ n)
  234.              (arraycall 'flonum (nsymbol-array '$besselarray) n))
  235.          (setf (arraycall 'flonum (nsymbol-array '$besselarray) k)
  236.                (simplify (list '(mplus)
  237.                        (simplify (list '(mtimes)
  238.                                '$%i
  239.                                (aref cyi k)))
  240.                        (aref cyr k)))))))))))
  241.  
  242. (defun bessel-k (arg order)
  243.   (cond ((zerop (imagpart arg))
  244.      ;; We have numeric args and the first arg is purely
  245.      ;; real. Call the real-valued Bessel function.  We call i0
  246.      ;; and i1 instead of jn, if possible.
  247.      (let ((arg (realpart arg)))
  248.        (cond ((zerop order)
  249.           (slatec:dbesk0 arg))
  250.          ((= order 1)
  251.           (slatec:dbesk1 arg))
  252.          (t
  253.           (multiple-value-bind (n alpha)
  254.               (floor (float order))
  255.             (let ((jvals (make-array (1+ n) :element-type 'double-float)))
  256.               (slatec:dbesk (float (realpart arg)) alpha 1 (1+ n) jvals 0)
  257.               (narray $besselarray $float n)
  258.               (fillarray (nsymbol-array '$besselarray) jvals)
  259.               (aref jvals n)))))))
  260.     (t
  261.      ;; The first arg is complex.  Use the complex-valued Bessel
  262.      ;; function
  263.      (multiple-value-bind (n alpha)
  264.          (floor (float order))
  265.        (let ((cyr (make-array (1+ n) :element-type 'double-float))
  266.          (cyi (make-array (1+ n) :element-type 'double-float)))
  267.          (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
  268.                     v-cyr v-cyi v-nz v-ierr)
  269.          (slatec::zbesk (float (realpart arg))
  270.                 (float (imagpart arg))
  271.                 alpha
  272.                 1
  273.                 (1+ n)
  274.                 cyr
  275.                 cyi
  276.                 0
  277.                 0)
  278.            (declare (ignore v-zr v-zi v-fnu v-kode v-n
  279.                 v-cyr v-cyi))
  280.  
  281.            ;; We should check for errors here based on the
  282.            ;; value of v-ierr.
  283.            (when (plusp v-ierr)
  284.          (format t "zbesk ierr = ~A~%" v-ierr))
  285.            (narray $besselarray $complete (1+ n))
  286.            (dotimes (k (1+ n)
  287.              (arraycall 'flonum (nsymbol-array '$besselarray) n))
  288.          (setf (arraycall 'flonum (nsymbol-array '$besselarray) k)
  289.                (simplify (list '(mplus)
  290.                        (simplify (list '(mtimes)
  291.                                '$%i
  292.                                (aref cyi k)))
  293.                        (aref cyr k)))))))))))
  294.  
  295. (defun $bessel_k (arg order)
  296.   (if (and (numberp order)
  297.        (numberp ($realpart arg))
  298.        (numberp ($imagpart arg)))
  299.       (bessel-k (complex ($realpart arg) ($imagpart arg)) order)
  300.       (subfunmakes '$bessel_k (ncons order) (ncons arg))))
  301.  
  302.  
  303.  
  304.  
  305. ;; I think g0(x) = exp(-x)*I[0](x), g1(x) = exp(-x)*I[1](x), and
  306. ;; gn(x,n) = exp(-x)*I[n](x), based on some simple numerical
  307. ;; evaluations.
  308.  
  309. (defun $g0 ($x)
  310.   (cond ((numberp $x)
  311.      (slatec:dbsi0e (float $x)))
  312.     (t (list '($g0 simp) $x))))
  313.  
  314. (defun $g1 ($x)
  315.   (cond ((numberp $x)
  316.      (slatec:dbsi1e (float $x)))
  317.     (t (list '($g1 simp) $x))))
  318.  
  319.  
  320. (declare-top (fixnum i n) (flonum x q1 q0 fn fi b1 b0 b an a1 a0 a)) 
  321.  
  322. (defun $gn ($x $n)
  323.   (cond ((and (numberp $x) (integerp $n))
  324.      (multiple-value-bind (n alpha)
  325.          (floor (float $n))
  326.        (let ((jvals (make-array (1+ n) :element-type 'double-float)))
  327.          (slatec:dbesi (float $x) alpha 2 (1+ n) jvals 0)
  328.          (narray $iarray $float n)
  329.          (fillarray (nsymbol-array '$iarray) jvals)
  330.          (aref jvals n))))
  331.     (t (list '(gn simp) $x $n))))
  332.  
  333.  
  334.  
  335. (declare-top(flonum rz cz a y $t t0 t1 d r1 rp sqrp rnpa r2 ta rn rl rnp rr cr rs cs rlam
  336.          clam qlam s phi rsum csum)
  337.      (fixnum n k1 k m mpo ln l ind)
  338.      (notype ($bessel notype notype) (bessel flonum flonum flonum))
  339.      (array* (flonum rj-bessel-array 1. cj-bessel-array 1.)
  340.          (notype $besselarray 1.))
  341.      (*fexpr $array))
  342.  
  343. ;; Bessel function of the first kind for real or complex arg and real
  344. ;; non-negative order.
  345. (defun $bessel ($arg $order)
  346.   (let ((a (float $order)))
  347.     (cond ((not (and (numberp $order)
  348.              (not (< a 0.0))
  349.              (numberp ($realpart $arg))
  350.              (numberp ($imagpart $arg))))
  351.        ;; Args aren't numeric.  Return unevaluated.
  352.        (list '($bessel simp) $arg $order))
  353.       ((zerop ($imagpart $arg))
  354.        ;; We have numeric args and the first arg is purely
  355.        ;; real. Call the real-valued Bessel function.  (Should we
  356.        ;; try calling j0 and j1 as appropriate instead of jn?)
  357.        (cond ((= $order 0)
  358.           (slatec:dbesj0 (float $arg)))
  359.          ((= $order 1)
  360.           (slatec:dbesj1 (float $arg)))
  361.          (t
  362.           (multiple-value-bind (n alpha)
  363.               (floor (float $order))
  364.             (let ((jvals (make-array (1+ n) :element-type 'double-float)))
  365.               ;; Use analytic continuation formula A&S 9.1.35:
  366.               ;;
  367.               ;; %j[v](z*exp(m*%pi*%i)) = exp(m*%pi*%i*v)*%j[v](z)
  368.               ;;
  369.               ;; for an integer m.  In particular, for m = 1:
  370.               ;;
  371.               ;; %j[v](-x) = exp(v*%pi*%i)*%j[v](x)
  372.               (cond ((>= $arg 0)
  373.                  (slatec:dbesj (float $arg) alpha (1+ n) jvals 0)
  374.                  (narray $besselarray $float n)
  375.                  (fillarray (nsymbol-array '$besselarray) jvals)
  376.                  (aref jvals n))
  377.                 (t
  378.                  (slatec:dbesj (- (float $arg)) alpha (1+ n) jvals 0)
  379.                  (narray $besselarray $complete n)
  380.                  (let ((s (cis (* $order pi))))
  381.                    (dotimes (k (1+ n))
  382.                  (let ((v (* s (aref jvals k))))
  383.                    (setf (arraycall 'flonum (nsymbol-array '$besselarray) k)
  384.                      (simplify `((mplus)
  385.                              ,(realpart v)
  386.                              ((mtimes)
  387.                               $%i
  388.                               ,(imagpart v)))))))
  389.                    (arraycall 'flonum (nsymbol-array '$besselarray) n)))))))))
  390.       (t
  391.        ;; The first arg is complex.  Use the complex-valued Bessel
  392.        ;; function
  393.        (multiple-value-bind (n alpha)
  394.            (floor (float $order))
  395.          (let ((cyr (make-array (1+ n) :element-type 'double-float))
  396.            (cyi (make-array (1+ n) :element-type 'double-float)))
  397.            (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
  398.                       v-cyr v-cyi v-nz v-ierr)
  399.            (slatec:zbesj (float ($realpart $arg))
  400.                  (float ($imagpart $arg))
  401.                  alpha
  402.                  1
  403.                  (1+ n)
  404.                  cyr
  405.                  cyi
  406.                  0
  407.                  0)
  408.          (declare (ignore v-zr v-zi v-fnu v-kode v-n
  409.                   v-cyr v-cyi v-nz))
  410.  
  411.          ;; Should check the return status in v-ierr of this
  412.          ;; routine.
  413.          (when (plusp v-ierr)
  414.            (format t "zbesj ierr = ~A~%" v-ierr))
  415.          (narray $besselarray $complete (1+ n))
  416.          (dotimes (k (1+ n)
  417.                (arraycall 'flonum (nsymbol-array '$besselarray) n))
  418.            (setf (arraycall 'flonum (nsymbol-array '$besselarray) k)
  419.              (simplify (list '(mplus)
  420.                      (simplify (list '(mtimes)
  421.                              '$%i
  422.                              (aref cyi k)))
  423.                      (aref cyr k))))))))))))
  424.  
  425. (defun $bessel_j (arg order)
  426.   (if (and (numberp order)
  427.        (numberp ($realpart arg))
  428.        (numberp ($imagpart arg)))
  429.       ($bessel (complex ($realpart arg) ($imagpart arg)) order)
  430.       (subfunmakes '$bessel_j (ncons order) (ncons arg))))
  431.  
  432. ;; Bessel function of the second kind, Y[n](z), for real or complex z
  433. ;; and non-negative real n.
  434. (defun bessel-y (arg order)
  435.   (cond ((zerop (imagpart arg))
  436.      ;; We have numeric args and the first arg is purely
  437.      ;; real. Call the real-valued Bessel function.
  438.      ;;
  439.      ;; For negative values, use the analytic continuation formula
  440.      ;; A&S 9.1.36:
  441.      ;;
  442.      ;; %y[v](z*exp(m*%pi*%i)) = exp(-v*m*%pi*%i) * %y[v](z)
  443.      ;;       + 2*%i*sin(m*v*%pi)*cot(v*%pi)*%j[v](z)
  444.      ;;
  445.      ;; In particular for m = 1:
  446.      ;;
  447.      ;; %y[v](-z) = exp(-v*%pi*%i) * %y[v](z) + 2*%i*cos(v*%pi)*%j[v](z)
  448.      ;; 
  449.      (let ((arg (realpart arg)))
  450.        (cond ((zerop order)
  451.           (cond ((>= arg 0)
  452.              (slatec:dbesy0 arg))
  453.             (t
  454.              ;; For v = 0, this simplifies to
  455.              ;;
  456.              ;; %y[0](-z) = %y[0](z) + 2*%i*%j[0](z)
  457.              (simplify `((mplus)
  458.                      ,(slatec:dbesy0 (- arg))
  459.                      ((mtimes)
  460.                       $%i
  461.                       ,(* 2 (slatec:dbesj0 (- arg)))))))))
  462.          ((= order 1)
  463.           (cond ((>= arg 0)
  464.              (slatec:dbesy1 arg))
  465.             (t
  466.              ;; For v = 1, this simplifies to
  467.              ;;
  468.              ;; %y[1](-z) = -%y[1](z) - 2*%i*%j[1](v)
  469.              (simplify `((mplus)
  470.                        ,(slatec:dbesy1 (- arg))
  471.                        ((mtimes)
  472.                     $%i
  473.                     ,(* -2 (slatec:dbesj1 (- arg)))))))))
  474.          (t
  475.           (multiple-value-bind (n alpha)
  476.               (floor (float order))
  477.             (let ((jvals (make-array (1+ n) :element-type 'double-float)))
  478.               (cond ((>= arg 0)
  479.                  (slatec:dbesy (float (realpart arg)) alpha (1+ n) jvals)
  480.                  (narray $besselarray $float n)
  481.                  (fillarray (nsymbol-array '$besselarray) jvals)
  482.                  (aref jvals n))
  483.                 (t
  484.                  (let* ((j ($bessel (- arg) order))
  485.                     (s1 (cis (- (* v pi))))
  486.                     (s2 (* #c(0 2) (cos (* v pi)))))
  487.                    (slatec:dbesy (- (float arg)) alpha (1+ n) jvals)
  488.                    (narray $yarray $complete n)
  489.                    (dotimes (k (1+ n))
  490.                  (let ((v (+ (* s1 (aref jvals k))
  491.                          (* s2 (arraycall 'flonum (nsymbol-array $besselarray) k)))))
  492.                    (setf (arraycall 'flonum (nsymbol-array '$yarray) k)
  493.                      (simplify `((mplus)
  494.                              ,(realpart v)
  495.                              ((mtimes)
  496.                               $%i
  497.                               ,(imagpart v)))))))
  498.                    (arraycall 'flonum (nsymbol-array '$yarray) n))))))))))
  499.     (t
  500.      ;; The first arg is complex.  Use the complex-valued Bessel
  501.      ;; function
  502.      (multiple-value-bind (n alpha)
  503.          (floor (float order))
  504.        (let ((cyr (make-array (1+ n) :element-type 'double-float))
  505.          (cyi (make-array (1+ n) :element-type 'double-float))
  506.          (cwrkr (make-array (1+ n) :element-type 'double-float))
  507.          (cwrki (make-array (1+ n) :element-type 'double-float)))
  508.          (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
  509.                     v-cyr v-cyi v-nz
  510.                     v-cwrkr v-cwrki v-ierr)
  511.          (slatec::zbesy (float (realpart arg))
  512.                 (float (imagpart arg))
  513.                 alpha
  514.                 1
  515.                 (1+ n)
  516.                 cyr
  517.                 cyi
  518.                 0
  519.                 cwrkr
  520.                 cwrki
  521.                 0)
  522.            (declare (ignore v-zr v-zi v-fnu v-kode v-n
  523.                 v-cyr v-cyi v-cwrkr v-cwrki))
  524.  
  525.            ;; We should check for errors here based on the
  526.            ;; value of v-ierr.
  527.            (when (plusp v-ierr)
  528.          (format t "zbesy ierr = ~A~%" v-ierr))
  529.            (narray $besselarray $complete (1+ n))
  530.            (dotimes (k (1+ n)
  531.              (arraycall 'flonum (nsymbol-array '$besselarray) n))
  532.          (setf (arraycall 'flonum (nsymbol-array '$besselarray) k)
  533.                (simplify (list '(mplus)
  534.                        (simplify (list '(mtimes)
  535.                                '$%i
  536.                                (aref cyi k)))
  537.                        (aref cyr k)))))))))))
  538.  
  539. (defun $bessel_y (arg order)
  540.   (if (and (numberp order)
  541.        (numberp ($realpart arg))
  542.        (numberp ($imagpart arg)))
  543.       (bessel-y (complex ($realpart arg) ($imagpart arg)) order)
  544.       (subfunmakes '$bessel_y (ncons order) (ncons arg))))
  545.  
  546. (declare-top(flonum rz y rs cs third sin60 term sum fi cossum sinsum sign (airy flonum)))
  547.  
  548. ;here is Ai'
  549. ;airy1(z):=if z = 0. then -1/(gamma(1/3.)*3.^(1/3.))
  550. ;else block([zz],z:-z,zz:2./3.*z^(3./2.),bessel(zz,4./3.),
  551. ;j:realpart(2/(3.*zz)*besselarray[0]-besselarray[1]),
  552. ;-1/3.*z*(j-realpart(bessel(zz,2./3.))));
  553.  
  554. (defun $airy ($arg)
  555.   (cond ((numberp $arg)
  556.      (slatec:dai (float $arg)))
  557.     (t
  558.      (list '($airy simp) $arg))))
  559.  
  560. (declare-top (flonum im re ys xs y x c t2 t1 s2 s1 s r2 r1 lamb h2 h)
  561.      (fixnum np1 n nu capn)
  562.      (notype (z-function flonum flonum))) 
  563.  
  564. (defun z-function (x y) 
  565.        ((lambda (xs ys capn nu np1 h h2 lamb r1 r2 s s1 s2 t1 t2 c bool re im) 
  566.         (setq xs (cond ((> 0.0 x) -1.0) (t 1.0)))
  567.         (setq ys (cond ((> 0.0 y) -1.0) (t 1.0)))
  568.         (setq x (abs x) y (abs y))
  569.         (cond ((and (> 4.29 y) (> 5.33 x))
  570.                (setq s (*$ (1+$ (*$ -0.23310023 y))
  571.                    (sqrt (1+$ (*$ -0.035198873 x x)))))
  572.                (setq h (*$ 1.6 s) h2 (*$ 2.0 h) capn (f+ 6. (fix (*$ 23.0 s))))
  573.                (setq nu (f+ 9. (fix (*$ 21.0 s)))))
  574.               (t (setq h 0.0) (setq capn 0.) (setq nu 8.)))
  575.         (and (> h 0.0) (setq lamb (^$ h2 capn)))
  576.         (setq bool (or (= h 0.0) (= lamb 0.0)))
  577.         (do ((n nu (f1- n)))
  578.             ((> 0. n))
  579.             (setq np1 (f1+ n))
  580.             (setq t1 (+$ h (*$ (float np1) r1) y))
  581.             (setq t2 (-$ x (*$ (float np1) r2)))
  582.             (setq c (//$ 0.5 (+$ (*$ t1 t1) (*$ t2 t2))))
  583.             (setq r1 (*$ c t1) r2 (*$ c t2))
  584.             (cond ((and (> h 0.0) (not (< capn n)))
  585.                (setq t1 (+$ s1 lamb) s1 (-$ (*$ r1 t1) (*$ r2 s2)))
  586.                (setq s2 (+$ (*$ r1 s2) (*$ r2 t1)) lamb (//$ lamb h2)))))
  587.         (setq im (cond ((= y 0.0) (*$ 1.77245384 (exp (-$ (*$ x x)))))
  588.                    (t (*$ 2.0 (cond (bool r1) (t s1))))))
  589.         (setq re (*$ -2.0 (cond (bool r2) (t s2))))
  590.         (cond ((> ys 0.0) (setq re (*$ re xs)))
  591.               (t (setq r1 (*$ 3.5449077 (exp (-$ (*$ y y) (*$ x x)))))
  592.              (setq r2 (*$ 2.0 x y))
  593.              (setq re (*$ (-$ re (*$ r1 (sin r2))) xs))
  594.              (setq im (-$ (*$ r1 (cos r2)) im))))
  595.         (list '(mlist simp) re im))
  596.     (cond ((> 0.0 x) -1.0) (t 1.0))
  597.     (cond ((> 0.0 x) -1.0) (t 1.0))
  598.     0. 0. 0. 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 nil 0.0 0.0)) 
  599.  
  600. (defun $nzeta ($z) 
  601.   (prog ($x $y $w) 
  602.     (cond ((and (numberp (setq $x ($realpart $z)))
  603.             (numberp (setq $y ($imagpart $z))))
  604.            (setq $w (z-function (float $x) (float $y)))
  605.            (return (simplify (list '(mplus)
  606.                        (simplify (list '(mtimes)
  607.                                (meval1 '$%i)
  608.                                (caddr $w)))
  609.                        (cadr $w)))))
  610.           (t (return (list '($nzeta simp) $z))))))
  611.  
  612.  
  613. (defun $nzetar ($z)
  614.   (prog ($x $y $w) 
  615.     (cond ((and (numberp (setq $x ($realpart $z)))
  616.             (numberp (setq $y ($imagpart $z))))
  617.            (setq $w (z-function (float $x) (float $y)))
  618.            (return (cadr $w)))
  619.           (t (return (list '($nzetar simp) $z))))))
  620.  
  621.  
  622. (defun $nzetai ($z)
  623.   (prog ($x $y $w) 
  624.     (cond ((and (numberp (setq $x ($realpart $z)))
  625.             (numberp (setq $y ($imagpart $z))))
  626.            (setq $w (z-function (float $x) (float $y)))
  627.            (return (caddr $w)))
  628.           (t (return (list '($nzetai simp) $z))))))
  629.  
  630.  
  631. ;; Initialize tables for Marsaglia's Ziggurat method of generating
  632. ;; random numbers.  See http://www.jstatsoft.org for a reference.
  633. ;;
  634. ;; Let 0 = x[0] < x[1] < x[2] <...< x[n].  Select a set of rectangles
  635. ;; with common area v such that
  636. ;;
  637. ;; x[k]*(f(x[k-1]) - f(x[k])) = v
  638. ;;
  639. ;; and
  640. ;;
  641. ;;              inf
  642. ;; v = r*f(r) + int f(x) dx
  643. ;;               r
  644. ;;
  645. ;; where r = x[n].
  646. ;;
  647. (defun ziggurat-init (n r v scale f finv)
  648.   ;; n = one less than the number of elements in the tables
  649.   ;; r = x[n]
  650.   ;; v = common area term
  651.   ;; scale = 2^scale is the scaling to use to make integers
  652.   ;; f = density function
  653.   ;; finv = inverse density function
  654.   (let ((x (make-array (1+ n) :element-type 'double-float))
  655.     (fx (make-array (1+ n) :element-type 'double-float))
  656.     (k-table (make-array (1+ n) :element-type '(unsigned-byte 32)))
  657.     (w-table (make-array (1+ n) :element-type 'double-float)))
  658.     (setf (aref x n) r)
  659.     (loop for k from (1- n) downto 1 do
  660.       (let ((prev (aref x (1+ k))))
  661.         (setf (aref x k) (funcall finv (+ (/ v prev)
  662.                           (funcall f prev))))
  663.         (setf (aref fx k) (funcall f (aref x k)))))
  664.  
  665.     (setf (aref x 0) 0d0)
  666.     (setf (aref fx 0) (funcall f (aref x 0)))
  667.     (setf (aref fx n) (funcall f (aref x n)))
  668.  
  669.     (loop for k from 1 to n do
  670.       (setf (aref k-table k)
  671.         (floor (scale-float (/ (aref x (1- k)) (aref x k)) scale)))
  672.       (setf (aref w-table k)
  673.         (* (aref x k) (expt .5d0 scale))))
  674.  
  675.     (setf (aref k-table 0) (floor (scale-float (/ (* r (funcall f r)) v) scale)))
  676.     (setf (aref w-table 0) (* (/ v (funcall f r)) (expt 0.5d0 scale)))
  677.     (values k-table w-table fx)))
  678.  
  679. ;; Marsaglia's Ziggurat method for Gaussians
  680. (let ((r 3.442619855899d0))
  681.   (flet ((density (x)
  682.        (declare (double-float x)
  683.             (optimize (speed 3) (safety 0)))
  684.        (exp (* -0.5d0 x x))))
  685.     (declaim (inline density))
  686.     (multiple-value-bind (k-table w-table f-table)
  687.     (ziggurat-init 127 r 9.91256303526217d-3 31
  688.                #'density
  689.                #'(lambda (x)
  690.                (sqrt (* -2 (log x)))))
  691.       (defun gen-gaussian-variate-ziggurat (state)
  692.     (declare (random-state state)
  693.          (optimize (speed 3)))
  694.     (loop
  695.         ;; We really want a signed 32-bit random number. So make a
  696.         ;; 32-bit unsigned number, take the low 31 bits as the
  697.         ;; number, and use the most significant bit as the sign.
  698.         ;; Doing this in other ways can cause consing.
  699.         (let* ((ran (random (ash 1 32) state))
  700.            (sign (ldb (byte 1 31) ran))
  701.            (j (if (plusp sign)
  702.               (- (ldb (byte 31 0) ran))
  703.               (ldb (byte 31 0) ran)))
  704.            (i (logand j 127))
  705.            (x (* j (aref w-table i))))
  706.           (when (< (abs j) (aref k-table i))
  707.         (return x))
  708.           (when (zerop i)
  709.         (loop
  710.             (let ((x (/ (- (log (random 1d0 state))) r))
  711.               (y (- (log (random 1d0 state)))))
  712.               (when (> (+ y y) (* x x))
  713.             (return-from gen-gaussian-variate-ziggurat
  714.               (if (plusp j)
  715.                   (- (+ r x))
  716.                   (+ r x)))))))
  717.           (when (< (* (random 1d0 state) (- (aref f-table (1- i))
  718.                         (aref f-table i)))
  719.                (- (density x) (aref f-table i)))
  720.         (return x))))))))
  721.  
  722. (defun $gauss ($mean $sd)
  723.   (cond ((and (numberp $mean) (numberp $sd))
  724.      (+$ (float $mean) (*$ (float $sd) (gen-gaussian-variate-ziggurat *random-state*))))
  725.     (t (list '($gauss simp) $mean $sd))))
  726.  
  727. (declare-top (flonum x w y (expint flonum)))
  728.  
  729. ;; I think this is the function E1(x).  At least some simple numerical
  730. ;; tests show that this expint matches the function de1 from SLATEC
  731.  
  732. ;; Exponential integral E1(x).  The Cauchy principal value is used for
  733. ;; negative x.
  734. (defun $expint (x)
  735.   (cond ((numberp x)
  736.      (values (slatec:de1 (float x))))
  737.     (t
  738.      (list '($expint simp) x))))
  739.  
  740.  
  741. ;; Define the Bessel funtion J[n](z)
  742.  
  743. (defprop $bessel_j bessel-j-simp specsimp)
  744.  
  745. (defprop $bessel_j
  746.     ((n x)
  747.      ((%derivative) ((mqapply) (($bessel_j array) n) x) n 1)
  748.      ((mplus)
  749.       ((mqapply) (($bessel_j array) ((mplus) -1 n)) x)
  750.       ((mtimes) -1 n ((mqapply) (($bessel_j array) n) x) ((mexpt) x -1))))
  751.   grad)
  752.  
  753. ;; If E is a maxima ratio with a denominator of DEN, return the ratio
  754. ;; as a Lisp rational.  Otherwise NIL.
  755. (defun max-numeric-ratio-p (e den)
  756.   (if (and (listp e)
  757.        (eq 'rat (caar e))
  758.        (= den (third e))
  759.        (integerp (second e)))
  760.       (/ (second e) (third e))
  761.       nil))
  762.  
  763. ;; For n = 1, computes f1(z) = 1/z*diff(f(z), z)
  764. ;; For n = 2, computes f2(z) = 1/z*diff(f1(z), z)
  765. ;; etc.
  766. (defun diffz (exp arg n)
  767.   (if (zerop n)
  768.       (simplify exp)
  769.       (diffz `((mtimes)
  770.            ((mexpt) ,arg -1)
  771.            ,(simplify ($diff exp arg)))
  772.          arg
  773.          (1- n))))
  774.  
  775. ;; Compute the Bessel function of half-integral order.
  776. ;;
  777. ;; ARG is the argument of the Bessel function
  778. ;; ORDER is the order, which must be half of an odd integer
  779. ;; 
  780. (defun bessel-jy-half-order (arg order pos-function neg-function)
  781.   ;; The description given here is for J functions, but they equally
  782.   ;; apply to Y.
  783.   ;;
  784.   ;; J[n+1/2](z) and J[-n-1/2](z) can be expressed in
  785.   ;; terms of elementary functions.  See A&S 9.1.30, let k = 1:
  786.   ;;
  787.   ;; diff(z^v * %j[v](z), z) = z^v * %j[v-1](z)
  788.   ;;
  789.   ;; and
  790.   ;;
  791.   ;; diff(z^(-v) * %j[v](z), z) = -z^(-v) * %j[v+1](z)
  792.   ;;
  793.   ;; We have
  794.   ;;
  795.   ;;   J[1/2](z) = sqrt(2/%pi)*sin(z)/sqrt(z)
  796.   ;;
  797.   ;; and
  798.   ;;   J[-1/2](z) = sqrt(2/%pi)*cos(z)/sqrt(z)
  799.   ;;
  800.   (cond ((plusp order)
  801.      ;; Setting v = 1/2 in the second formula of A&S 9.1.30 we have
  802.      ;;
  803.      ;; J[n+1/2](z) = (-1)^n * sqrt(z) * z^n *
  804.      ;;                diffz(J[1/2](z),z,n)
  805.      ;;
  806.      ;; where diffz(f,z,n) means compute 1/z*diff(f,z) n times.
  807.      ;;
  808.      ;; or, using the expressin for J[1/2](z) above:
  809.      ;;
  810.      (let* ((n (floor order))
  811.         (var (gensym))
  812.         ;; deriv = diff(sin(z)/z,z,n)
  813.         (deriv (subst arg
  814.                   var
  815.                   (diffz `((mtimes simp)
  816.                        ((mexpt simp) ,var -1)
  817.                        ((,pos-function simp) ,var))
  818.                      var
  819.                      n))))
  820.        (simplify `((mtimes)
  821.                ((mexpt) 2 ((rat) 1 2))
  822.                ((mexpt) $%pi ((rat) -1 2))
  823.                ((mexpt) -1 ,n)
  824.                ((mexpt) ,arg ,n)
  825.                ((mexpt) ,arg ((rat) 1 2))
  826.                ,deriv))))
  827.     (t
  828.      ;; We use the first formula and J[-1/2](z) above to get
  829.      ;;
  830.      ;; J[-n-1/2](z) = sqrt(z) * sqrt(2/%pi) * z^n
  831.      ;;                    diffz(cos(z)/z, z, n);
  832.      (let* ((n (floor (- order)))
  833.         (var (gensym))
  834.         ;; deriv = diff(cos(z)/z,z,n)
  835.         (deriv (subst arg var
  836.                   (diffz `((mtimes simp)
  837.                        ((mexpt simp) ,var -1)
  838.                        ((,neg-function simp) ,var))
  839.                      var
  840.                      n))))
  841.        (simplify `((mtimes)
  842.                ((mexpt) 2 ((rat) 1 2))
  843.                ((mexpt) $%pi ((rat) -1 2))
  844.                ((mexpt) ,arg ,n)
  845.                ((mexpt) ,arg ((rat) 1 2))
  846.                ,deriv))))))
  847.  
  848. (defun bessel-i-half-order (arg order pos-function neg-function)
  849.   ;; I[n+1/2](z) and I[-n-1/2](z) can be expressed in
  850.   ;; terms of elementary functions.  See A&S 9.6.28:
  851.   ;;
  852.   ;;
  853.   ;;          k
  854.   ;; [ 1   d ]
  855.   ;; [--- ---]  [z^v*I[v](z)] = z^(v-k) * I[v-k](z)
  856.   ;; [ z   dz]
  857.   ;;
  858.   ;;
  859.   ;;          k
  860.   ;; [ 1   d ]
  861.   ;; [--- ---]  [z^(-v)*I[v](z)] = z^(-v-k) * I[v+k](z)
  862.   ;; [ z   dz]
  863.   ;;
  864.   ;;
  865.   ;; We have
  866.   ;;
  867.   ;;   I[1/2](z) = sqrt(2/%pi)*sinh(z)/sqrt(z)
  868.   ;;
  869.   ;;   I[-1/2](z) = sqrt(2/%pi)*cosh(z)/sqrt(z)
  870.   ;;
  871.   ;; Thus, for the second equation above, v = 1/2
  872.   ;;
  873.   ;;          k
  874.   ;; [ 1   d ]
  875.   ;; [--- ---]  [1/sqrt(z)* sqrt(2/%pi) * sinh(z)/sqrt(z)] = 1/sqrt(z)/z^k * I[1/2+k](z)
  876.   ;; [ z   dz]
  877.   ;;
  878.   ;; or
  879.   ;;
  880.   ;;          k
  881.   ;; [ 1   d ]
  882.   ;; [--- ---]  [sqrt(2/%pi)*sinh(z)/z] = 1/sqrt(z)/z^k * I[1/2+k](z)
  883.   ;; [ z   dz] 
  884.   ;;
  885.   ;;
  886.   ;; So
  887.   ;;                                                    k
  888.   ;;                                           [ 1   d ]
  889.   ;; I[1/2+k](z) = sqrt(z) * z^k * sqrt(2/%pi) [--- ---]  [sinh(z)/z]
  890.   ;;                                           [ z   dz]
  891.   ;;
  892.   ;;
  893.   ;; Similarly, for the first equation above with v = -1/2
  894.   ;;
  895.   ;;          k
  896.   ;; [ 1   d ]
  897.   ;; [--- ---]  [1/sqrt(z)* sqrt(2/%pi) * cosh(z)/sqrt(z)] = 1/sqrt(z)/z^k * I[-1/2-k](z)
  898.   ;; [ z   dz]
  899.   ;;
  900.   ;; or
  901.   ;;                                                    k
  902.   ;;                                            [ 1   d ]
  903.   ;; I[-1/2-k](z) = sqrt(z) * z^k * sqrt(2/%pi) [--- ---]  [cosh(z)/z]
  904.   ;;                                            [ z   dz]
  905.  
  906.   (let* ((n (truncate (abs order)))
  907.      (var (gensym))
  908.      (deriv (subst arg var
  909.                (simplify (diffz `((mtimes simp)
  910.                       ((mexpt simp) ,var -1)
  911.                       ((,(if (plusp order) pos-function neg-function)
  912.                         simp) ,var))
  913.                     var
  914.                     n)))))
  915.     (simplify `((mtimes)
  916.         ((mexpt) 2 ((rat) 1 2))
  917.         ((mexpt) $%pi ((rat) -1 2))
  918.         ((mexpt) ,arg ((rat) 1 2))
  919.         ((mexpt) ,arg ,n)
  920.         ,deriv))))
  921.  
  922. (defun bessel-j-simp (exp ignored z)
  923.   (declare (ignore ignored))
  924.   (let ((order (simpcheck (car (subfunsubs exp)) z))
  925.     (rat-order nil))
  926.     (subargcheck exp 1 1 '$bessel_j)
  927.     (let* ((arg (simpcheck (car (subfunargs exp)) z)))
  928.       (when (and (numberp arg) (zerop arg)
  929.          (numberp order))
  930.     ;; J[v](0) = 1 if v = 0.  Otherwise 0.
  931.     (return-from bessel-j-simp
  932.       (if (zerop order)
  933.           1
  934.           0)))
  935.       (cond ((and $numer (numberp order)
  936.           (complex-number-p arg))
  937.          ;; We have numeric order and arg, so let's try to
  938.          ;; evaluate it numerically.
  939.          (let ((real-arg ($realpart arg))
  940.            (imag-arg ($imagpart arg)))
  941.            (cond ((or (and (floatp real-arg) (numberp imag-arg))
  942.               (and $numer (numberp real-arg) (numberp imag-arg)))
  943.               ;; Numerically evaluate it if the arg is a float
  944.               ;; or if the arg is a number and $numer is
  945.               ;; non-NIL.
  946.               ($bessel arg order))
  947.              ((minusp order)
  948.               ;; A&S 9.1.5
  949.               ;; J[-n](x) = (-1)^n*J[n](x)
  950.               (if (evenp order)
  951.               (subfunmakes '$bessel_j (ncons (- order)) (ncons arg))
  952.               `((mtimes simp) -1 ,(subfunmakes '$bessel_j (ncons (- order)) (ncons arg)))))
  953.              (t
  954.               (eqtest (subfunmakes '$bessel_j (ncons order) (ncons arg))
  955.                   exp)))))
  956.         ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
  957.          ;; When order is a fraction with a denominator of 2, we
  958.          ;; can express the result in terms of elementary
  959.          ;; functions.
  960.          ;;
  961.          ;; J[1/2](z) is a function of sin.  J[-1/2](z) is a
  962.          ;; function of cos.
  963.          (bessel-jy-half-order arg rat-order '%sin '%cos))
  964.         (t
  965.          (eqtest (subfunmakes '$bessel_j (ncons order) (ncons arg))
  966.              exp))))))
  967.  
  968.  
  969. ;; Define the Bessel funtion Y[n](z)
  970.  
  971. (defprop $bessel_y bessel-y-simp specsimp)
  972.  
  973. (defprop $bessel_y
  974.     ((n x)
  975.      ((%derivative) ((mqapply) (($bessel_y array) n) x) n 1)
  976.      ((mplus)
  977.       ((mqapply) (($bessel_y array) ((mplus) -1 n)) x)
  978.       ((mtimes) -1 n ((mqapply) (($bessel_y array) n) x) ((mexpt) x -1))))
  979.   grad)
  980.  
  981. (defun bessel-y-simp (exp ignored z)
  982.   (declare (ignore ignored))
  983.   (let ((order (simpcheck (car (subfunsubs exp)) z))
  984.     (rat-order nil))
  985.     (subargcheck exp 1 1 '$bessel_y)
  986.     (let* ((arg (simpcheck (car (subfunargs exp)) z)))
  987.       (cond ((and $numer (numberp order)
  988.           (complex-number-p arg))
  989.          (let ((real-arg ($realpart arg))
  990.            (imag-arg ($imagpart arg)))
  991.            (cond ((or (and (floatp real-arg) (numberp imag-arg))
  992.               (and $numer (numberp real-arg) (numberp imag-arg)))
  993.               ;; Numerically evaluate it if the arg is a float
  994.               ;; or if the arg is a number and $numer is
  995.               ;; non-NIL.
  996.               (bessel-y (complex real-arg imag-arg) (float order)))
  997.              ((minusp order)
  998.               ;; A&S 9.1.5
  999.               ;; Y[-n](x) = (-1)^n*Y[n](x)
  1000.               (if (evenp order)
  1001.               (subfunmakes '$bessel_y (ncons (- order)) (ncons arg))
  1002.               `((mtimes simp) -1 ,(subfunmakes '$bessel_y (ncons (- order)) (ncons arg)))))
  1003.              (t
  1004.               (eqtest (subfunmakes '$bessel_y (ncons order) (ncons arg))
  1005.                   exp)))))
  1006.         ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
  1007.          ;; When order is a fraction with a denominator of 2, we
  1008.          ;; can express the result in terms of elementary
  1009.          ;; functions.
  1010.          ;;
  1011.          ;; Y[1/2](z) = -J[1/2](z) is a function of sin.
  1012.          ;; Y[-1/2](z) = -J[-1/2](z) is a function of cos.
  1013.          (simplify `((mtimes) -1 ,(bessel-jy-half-order arg rat-order '%sin '%cos))))
  1014.         (t
  1015.          (eqtest (subfunmakes '$bessel_y (ncons order) (ncons arg))
  1016.              exp))))))
  1017.  
  1018. ;; Define the Bessel funtion I[n](z)
  1019.  
  1020. (defprop $bessel_i bessel-i-simp specsimp)
  1021.  
  1022. (defprop $bessel_i
  1023.     ((n x)
  1024.      ((%derivative) ((mqapply) (($bessel_i array) n) x) n 1)
  1025.      ((mplus)
  1026.       ((mqapply) (($bessel_i array) ((mplus) -1 n)) x)
  1027.       ((mtimes) -1 n ((mqapply) (($bessel_i array) n) x) ((mexpt) x -1))))
  1028.   grad)
  1029.  
  1030. (defun bessel-i-simp (exp ignored z)
  1031.   (declare (ignore ignored))
  1032.   (let ((order (simpcheck (car (subfunsubs exp)) z))
  1033.     (rat-order nil))
  1034.     (subargcheck exp 1 1 '$bessel_i)
  1035.     (let* ((arg (simpcheck (car (subfunargs exp)) z)))
  1036.       (cond ((and $numer (numberp order)
  1037.           (complex-number-p arg))
  1038.          (let ((real-arg ($realpart arg))
  1039.            (imag-arg ($imagpart arg)))
  1040.            (cond ((or (and (floatp real-arg) (numberp imag-arg))
  1041.               (and $numer (numberp real-arg) (numberp imag-arg)))
  1042.               ;; Numerically evaluate it if the arg is a float
  1043.               ;; or if the arg is a number and $numer is
  1044.               ;; non-NIL.
  1045.               (bessel-i (complex real-arg imag-arg) (float order)))
  1046.              ((minusp order)
  1047.               ;; A&S 9.6.6
  1048.               ;; I[-n](x) = I[n](x)
  1049.               (subfunmakes '$bessel_i (ncons (- order)) (ncons arg)))
  1050.              (t
  1051.               (eqtest (subfunmakes '$bessel_i (ncons order) (ncons arg))
  1052.                   exp)))))
  1053.         ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
  1054.          ;; When order is a fraction with a denominator of 2, we
  1055.          ;; can express the result in terms of elementary
  1056.          ;; functions.
  1057.          ;;
  1058.          ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z)
  1059.          ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z)
  1060.          (bessel-i-half-order arg rat-order '%sinh '%cosh))
  1061.         (t
  1062.          (eqtest (subfunmakes '$bessel_i (ncons order) (ncons arg))
  1063.              exp))))))
  1064.  
  1065. ;; Define the Bessel funtion K[n](z)
  1066.  
  1067. (defprop $bessel_k bessel-k-simp specsimp)
  1068.  
  1069. (defprop $bessel_k
  1070.     ((n x)
  1071.      ((%derivative) ((mqapply) (($bessel_k array) n) x) n 1)
  1072.      ((mplus simp)
  1073.       ((mtimes simp) -1 ((mqapply simp) (($bessel_k simp array) n) x))
  1074.       ((mtimes simp) -1 n ((mexpt simp) x -1)
  1075.        ((mqapply simp) (($bessel_k simp array) n) x))))
  1076.   grad)
  1077.  
  1078. (defun bessel-k-half-order (arg order)
  1079.   ;; K[n+1/2](z) and K[-n-1/2](z) can be expressed in terms of
  1080.   ;; elementary functions.  See A&S 9.6.28.  Let G[v](z) =
  1081.   ;; exp(v*%pi*%i)*K[v](z).
  1082.   ;;
  1083.   ;;          k
  1084.   ;; [ 1   d ]
  1085.   ;; [--- ---]  [z^v*G[v](z)] = z^(v-k) * G[v-k](z)
  1086.   ;; [ z   dz]
  1087.   ;;
  1088.   ;;
  1089.   ;;          k
  1090.   ;; [ 1   d ]
  1091.   ;; [--- ---]  [z^(-v)*G[v](z)] = z^(-v-k) * G[v+k](z)
  1092.   ;; [ z   dz]
  1093.   ;;
  1094.   ;;
  1095.   ;; or
  1096.   ;;          k
  1097.   ;; [ 1   d ]
  1098.   ;; [--- ---]  [z^v*K[v](z)] = z^(v-k) * K[v-k](z) * exp(k*%pi*%i)
  1099.   ;; [ z   dz]
  1100.   ;;
  1101.   ;;
  1102.   ;;          k
  1103.   ;; [ 1   d ]
  1104.   ;; [--- ---]  [z^(-v)*K[v](z)] = z^(-v-k) * K[v+k](z) * exp(k*%pi*%i)
  1105.   ;; [ z   dz]
  1106.   ;;
  1107.   ;;
  1108.   ;;
  1109.   ;; We have
  1110.   ;;
  1111.   ;;   K[1/2](z) = sqrt(2/%pi/z)*exp(-z) = K[-1/2](z)
  1112.   ;;
  1113.   ;; From the second equation, we have, with v = 1/2:
  1114.   ;;
  1115.   ;;          k
  1116.   ;; [ 1   d ]
  1117.   ;; [--- ---]  [1/sqrt(z) * sqrt(2/%pi/z) * exp(-z)] = 1/sqrt(z)/z^k * K[1/2+k](z) * exp(k*%pi*%i)
  1118.   ;; [ z   dz]
  1119.   ;; 
  1120.   ;; or
  1121.   ;;
  1122.   ;;                        k
  1123.   ;;               [ 1   d ]
  1124.   ;; K[1/2+k](z) = [--- ---]  [exp(-z)/z] * sqrt(2/%pi) * sqrt(z) * z^k * (-1)^k
  1125.   ;;               [ z   dz]
  1126.  
  1127.   
  1128.   (let* ((n (floor (abs order)))
  1129.      (var (gensym))
  1130.      ;; deriv = diff(exp(-z)/z,z,n)
  1131.      (deriv (subst arg var
  1132.                ($diff `((mtimes simp)
  1133.                 ((mexpt simp) ,var -1)
  1134.                 ((mexpt simp) $%e ((mtimes simp) -1 ,var)))
  1135.                   var
  1136.                   n))))
  1137.     (simplify `((mtimes)
  1138.         ,(if (evenp n) 1 -1)
  1139.         ((mexpt) 2 ((rat) 1 2))
  1140.         ((mexpt) $%pi ((rat) -1 2))
  1141.         ((mexpt) ,arg ((rat) 1 2))
  1142.         ((mexpt) ,arg ,n)
  1143.         ,deriv))))
  1144.   
  1145. (defun bessel-k-simp (exp ignored z)
  1146.   (declare (ignore ignored))
  1147.   (let ((order (simpcheck (car (subfunsubs exp)) z))
  1148.     (rat-order nil))
  1149.     (subargcheck exp 1 1 '$bessel_k)
  1150.     (let* ((arg (simpcheck (car (subfunargs exp)) z)))
  1151.       (cond ((and $numer (numberp order)
  1152.           (complex-number-p arg))
  1153.          (let ((real-arg ($realpart arg))
  1154.            (imag-arg ($imagpart arg)))
  1155.            (cond ((or (and (floatp real-arg) (numberp imag-arg))
  1156.               (and $numer (numberp real-arg) (numberp imag-arg)))
  1157.               ;; Numerically evaluate it if the arg is a float
  1158.               ;; or if the arg is a number and $numer is
  1159.               ;; non-NIL.
  1160.               (bessel-k (complex real-arg imag-arg) (float order)))
  1161.              ((minusp order)
  1162.               ;; A&S 9.6.6
  1163.               ;; K[-v](x) = K[v](x)
  1164.               (subfunmakes '$bessel_k (ncons (- order)) (ncons arg)))
  1165.              (t
  1166.               (eqtest (subfunmakes '$bessel_k (ncons order) (ncons arg))
  1167.                   exp)))))
  1168.         ((and $besselexpand
  1169.           (setq rat-order (max-numeric-ratio-p order 2)))
  1170.          ;; When order is a fraction with a denominator of 2, we
  1171.          ;; can express the result in terms of elementary
  1172.          ;; functions.
  1173.          ;;
  1174.          ;; K[1/2](z) = sqrt(2/%pi/z)*exp(-z) = K[1/2](z)
  1175.          (bessel-k-half-order arg rat-order))
  1176.         (t
  1177.          (eqtest (subfunmakes '$bessel_k (ncons order) (ncons arg))
  1178.              exp))))))
  1179.